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Hydrodynamical calculations have been successful in describing global observables in ultrarela- 
tivistic heavy ion collisions, which aim to observe the production of the quark-gluon plasma. On the 
other hand, recently, a lot of evidence that there exists a critical end point (CEP) in the QCD phase 
diagram has been accumulating. Nevertheless, so far, no equation of state with the CEP has been 
employed in hydrodynamical calculations. In this paper, we construct the equation of state with 
the CEP on the basis of the universality hypothesis and show that the CEP acts as an attractor of 
isentropic trajectories. We also consider the time evolution in the case with the CEP and discuss 
how the CEP affects the final state observables, such as the correlation length, fluctuation, chemical 
freezeout, kinetic freezeout, and so on. Finally, we argue that the anomalously low kinetic freezeout 
temperature at the BNL Relativistic Heavy Ion Collider suggests the possibility of the existence of 
the CEP. 



I. INTRODUCTION 

The structure of the QCD phase diagram is one of 
the most interesting and important topics in nuclear and 
particle physics, and a lot of intensive studies have been 
carried out from both theoretical side and experimen- 
tal side 0. Recent theoretical studies have revealed the 
rich structure in the QCD phase diagram at finite den- 
sity, which consists of not only the quark-gluon plasma 
(QGP) phase and hadron phase but also the two flavor 
color superconductor phase and color-flavor locked phase 
■ Furthermore, in the finite temperature direction, re- 
cent lattice calculations suggest the existence of mesonic 
bound states even above the deconfinement phase tran- 
sition temperature H S S H S Q . 

In this paper, we focus on the critical end point (CEP) 
in the QCD phase diagram, which is the terminating 
point of the first order phase transition . In partic- 
ular, we will consider its effects on the observables in 
ultrarelativistic heavy ion collisions in detail. Several ef- 
fective theory analyses have predicted the existence of 
the CEP in the QCD phase diagram [HI El El- Since 
Fodor and Katz presented a QCD phase diagram with 
the CEP for the first time with lattice calculation [13|, a 
lot of remarkable progress in finite temperature and den- 
sity lattice calculation has been made |14Lll5|. However, 
unfortunately, it will need more time before we reach the 
conclusive result on the existence of the CEP, its precise 
position, and so on, because finite density lattice QCD 
still has difficulties in performing accurate calculation on 
large lattices and at large chemical potential. Here, the 



1 The terminology 'critical end point' is sometimes used for a dif- 
ferent meaning [jj. However, it is rather customary in hadron 
physics to use it for the meaning we employ here. Thus, we use 
'critical end point' instead of 'critical point' in this paper. We 
thank M. Stephanov for informing us of Ref. 0. 



important point is that the existence of the CEP was 
shown semi-quantitatively by non-perturbative QCD cal- 
culation. The precise determination of the location of the 
CEP in the QCD phase diagram is a difficult problem. 
The locations of the CEP in different studies are, in fact, 
scattered over the T (temperature)-/^ (baryon chemical 
potential) plane 12]. Therefore, analyses from the exper- 
imental side are also indispensable in order to understand 
the properties of the QCD critical end point. 

A promising way to investigate the existence and lo- 
cation of the QCD critical end point is to carry out ex- 
periments at various collision energies and compare their 
results. However, we have to find the observables which 
are suitable for detecting the CEP [H El El El HJ 
for that purpose. In this paper, we study the CEP from 
the point of view of phenomenological analyses which can 
be compared with experimental results easily and discuss 
the indication of the existence of the CEP in experimen- 
tal observables. 

We divide the discussion into three steps. First, we 
construct equations of state that include the CEP. This 
step is necessary in order to consider the hydrodynami- 
cal evolution of the system. The guiding principle here is 
the theory of critical phenomena, in particular, the uni- 
versality hypothesis. Second, we consider the deviation 
from the thermal equilibrium near the CEP. This is im- 
portant because i) the typical time scale is elongated near 
the CEP, and ii) the system produced in ultrarelativistic 
heavy ion collisions cannot stay near the CEP due to the 
strong expansion. Finally, we discuss the possibility to 
observe the consequences of the CEP experimentally by 
using the preceding results. 

The structure of this paper is as follows. In Section 
II we argue the general feature of the equation of state 
with the CEP on the basis of the universality hypothesis, 
the focusing effect in isentropic trajectories in the T-fiB 
plane, and hydrodynamical expansion near the CEP. In 
Section III we investigate the growth and decrease of the 
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correlation length near the CEP from the point of view of 
slowing-out-of-equilibrium. In Section IV we discuss the 
consequences of the CEP in ultrarelativistic heavy ion 
collision experiments. Section V is devoted to summary. 



II. EQUATION OF STATE WITH THE 
CRITICAL END POINT 

In this section, we construct realistic equations of state 
including the CEP. The equation of state with the CEP 
consists of two parts, the singular part and non-singular 
part. We assume that the CEP in QCD belongs to the 
same universality class as that in the three dimensional 
Ising model on the basis of the universality hypothesis. 
After mapping the variables and the equation of state 
near the CEP in the three dimensional Ising model onto 
those in QCD, we match the singular entropy density 
near the CEP with the non-singular QGP and hadron 
phase entropy densities which are known away from the 
CEP. From this procedure we determine the behavior of 
the entropy density which includes both the singular part 
and non-singular part in a large region in the T-^lb plane. 
Finally, from the entropy density, we extract the behavior 
of other thermodynamical quantities such as the baryon 
number density, pressure, and energy density, and show 
that in the T-\ib plane the CEP acts as an attractor 
of isentropic trajectories, Hb/s = const., with ns and 
s being the baryon number density and entropy density, 
respectively. 



A. Singular Part of Equation of State 

In the three dimensional Ising model, the magnetiza- 
tion M (the order parameter) is a function of the reduced 
temperature r = (T — T c )/T c and the external magnetic 
field h with T c being the critical temperature. The CEP 
is located at the origin (r,h) = (0,0). At r < the 
order of the phase transition is first and at r > it is 
crossover. A useful form of the equation of state of the 
three-dimensional Ising model is given through the para- 
metric representation by the two variables R and 9, 



G(h,r) 



M = M O R 6, 
h = h o R 0S h(9) 

= h Q R ps {6 - 0.7620K9 3 
r = R(l-9 2 ) (R>0,- 



-0.008046* 5 ), 
•1.154 < 9 < 1.154), 



(1) 



where Mq and ho are normalization constants and crit- 
ical exponents (3 and S are 0.326 and 4.80, respectively 
|2lj . We set the normalization constants Mo and ho by 
imposing M(r=—l,h=+0) = 1 and M(r=0,h=l) = 1, 
which assures that the critical magnetization behaves as 
M(r = 0,h) cx sgn(/i)|/i| 1/<5 and M(r,h = +0) oc \rf 
(r < 0) around the origin. 

In order to determine the singular part of the en- 
tropy density, we start from the Gibbs free energy density 



G(h, r) = F(M, r) - Mh, 



(2) 



where F(M, r) is the free energy density. Setting the free 
energy density, 



F(M,r)=hoM o r 2 - a g(0), 



(3) 



together with the relation h = (dF/dM) r , we obtain the 
differential equation for g(6), 

h{9){\-9 2 + 2(38 2 ) = 2(2 - a)9 g (9) + (l-9 2 )g' (9), (4) 

where a critical exponent a is 0.11 |2l|. Substituting h(9) 
into Eq. we obtain 

g(9) = g(l) + c 1 (l-9 2 ) + c 2 (l - 9 2 ) 2 +c 3 (l ~9 2 ) 3 , (5) 

where g(l) is determined by Eqs. and Q, 
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Mi), 



2- a 

and the coefficients c±, c%, and C3 are given by 
1 1 



(6) 



Cl 



2 a - 1 
1 



{(1 - 2/3)(l + a + b) - 2(3{a + 2b)} , 



c 2 = - —{2/36- (1- 2/3) (a + 26)}. 
la 



1 



C3 



2(a+l) 



6(1-2/3), 



(J) 



where a = —0.76201 and 6 = 0.00804 are the coefficients 
of the 9 3 and 9 5 terms in h(9), respectively [2l|. Actually, 
there is freedom to add the function C(l — 9 2 ) 2 ~ a to 
g(9), with C being a constant. However, we will not take 
into account this part with a free parameter, because, 
whereas the universality tells us about only the critical 
behavior near the phase transition, this part does not 
show a singular behavior. 

Differentiating the Gibbs free energy by the tempera- 
ture, we obtain the singular part of the entropy density 
s c near the QCD critical end point, 



dG 
df 

dG\ dh /dG\ dr 



dhJ„dT \drj h dT' 



(8) 



where (757^) and (^r) h are, respectively, given by 



dG\ fdF(M,r)\ fdM\ , w 

Had*-" 

= ~M, (9) 
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,9rJ h V * Jh \ dr Jh 

= . h ° M ° R " _ x f(2-a)h'(0)g(9) - BOh' (0)1(6) - 86h(9)g'(9) + p5h 2 (9)\. (10) 

2859h(6) + (1 - 9 2 )h'(9) I i 




FIG. 1: Sketch of the r-h axes (three dimensional Ising 
model) mapped onto the T-/j,b plane (QCD). The r axis is 
tangential to the QCD phase boundary at the CEP. We set 
the h direction perpendicular to the r direction. 




FIG. 2: Dimensionless variable S c (T,fj,B) as a func- 
tion of T and [ib- The CEP is located at (T, (ab) = 
(154.7 MeV, 367.8 MeV) (the triangle in the contour plot). 



B. Thermodynamical Quantities 



Note that T in Eq. © is the temperature on the QCD 
side. In Eqs. © and H10|) . we use the differential equation 
Eq. 10} and the following relations, 



dh), 



dh / 
dR s 
dr /h 

dr , 



1 
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hoR? 6 - 1 2B59h(9) + (1 - 9 2 )h'(9) ' 

_L 1~9 2 

hvR? 5 2359h(8) + (1 - 9 2 )h'(9) ' 

h'(0) 

2/350h(0) + (l-6 2 )h'(6)' 

_{36 h(9) 

~~R 2359h(6) + (l-6 2 )h'(9)' 



(11) 



For quantitative construction of equations of state with 
the CEP, we fix the relation between the scales in (r, h) 
and (T, /ib) variables, which provides the size of the crit- 
ical region around the CEP in the T-fig plane, as follows: 
Ar = 1 (r-h plane) «-» A/is crit (T-^b plane) and Ah = 1 
(r-h plane) <-> AT crit (T-fj, B plane) (Fig.[TI). A^ Bcrit and 
AT cr it give, respectively, approximate extensions of the 
critical region in the \i b and T directions when the loca- 
tion of the CEP is close to the T axis and, as a result, 
the r axis is approximately parallel to the ub axis, as re- 
cent lattice calculations suggest [lj, [Tj, [l£j • In order to 
connect the equations of state in the singular region and 
the non-singular region smoothly, we define the dimen- 
sionless variable S c (T,[Ib) for the singular part of the 
entropy density s c , which has the dimension [energy] -1 , 



S c (T,n B ) = A (AT crit , A^ Bcrit )s c (T, n B ), 



(12) 



Now we map the r-h plane in the three dimensional 
Ising model onto the T-fig plane in QCD in order to 
determine dh/dT and dr/dT in Eq. ©. The CEP in 
the three dimensional Ising Model, which is the origin in 
the r-h plane, is mapped to the CEP in QCD, (T, /j,b) — 
(Te, /j-be)- The r axis is tangential to the first order 
phase transition line at the CEP (23. However, there 
is no general rule about how the h axis is mapped in 
the T-fiB plane. Here, for simplicity, we set the h axis 
perpendicular to the r axis. In Fig.^ at r < the order 
of the phase transition is first and at r > it is crossover. 



where A(AT critl A/x Bcrit ) = y/AT 2 ^ + A/i| crit x D and 
D is a dimensionless constant. The extension of the crit- 
ical domain around the CEP is specified by the param- 
eters AT cr ;t, Afj,Bcriti and D. In Fig. [3 an example of 
S c is shown as a function of T and /is. The CEP is 
located at (T,fi B ) =(157.4 MeV, 367.8 MeV) (the trian- 
gle in the contour plot). We use the same parameters 
(AT cr it, A/iB cr it, D) as those in Fig. ^ (left). S c is ob- 
tained by solving numerically the implicit relations de- 
rived in the preceding subsection. It is clearly seen that 
the order of the phase transition changes from crossover 
to first order at the CEP. The contour lines, S c = ±0.5, 
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give a rough idea about the size of the phase transition 
region. At lower /ig below the critical chemical poten- 
tial, the phase transition is smeared and the effect of 
the phase transition is observed in a larger domain in 



the temperature direction. Using the dimensionless vari- 
able S C (T, /is), we define the entropy density in the T-fig 
plane. 



s(T,Hb) = 7: (1 



tanh[S" c (r, fi B )]) sn(T,fj, B ) + - (1 



tanh^ (T, \i B )] ) s Q (T, \i B ) , 



(13) 



r 



where sh and sq are the entropy densities in the hadron 
phase and QGP phase away from the CEP, respectively. 
This entropy density includes both singular and non- 
singular contributions, and more importantly, gives the 
correct critical exponents near the QCD critical end 
point, sh is calculated from the equation of state of the 
hadron phase in the excluded volume approximation |23j| , 

i 

= ]>> idcal (T,^), (14) 



In the first order phase transition region, in order to 
take into account the discontinuity in the entropy den- 
sity and baryon number density on the phase boundary, 
it is necessary to add the following term to Eq. I|16l) for 

T > T c (fi B ), 



dT c ( PB ) 



dm 



(s(T c (mb)+0,/x b )-s(T c ( Ms )-0,Mb)), (17) 



where T c (/i B ) is the temperature of the first order phase 
transition at /i B . The pressure P is obtained by integrat- 
ing the entropy density with regard to the temperature, 



where P is the pressure, p^ doal is the pressure of the ideal 
gas of particle i, fi B i is the baryon chemical potential of 
particle i, and Vq is the hadron volume common to all 
hadron species. Non-strange resonances with mass up to 
2 GeV are included in the sum and the radius of hadrons 
is fixed at 0.7 fm. We obtain sq from the equation of 
state of the QGP phase in the Bag model, 



(32 + 2EV / V t4 



+ 



180 
4tt 2 V 3 J 



2 V 3 ) 



B, 



rj-i'2 



(15) 



where the number of the flavors Nf is 2 and the bag con- 
stant B is (220 MeV) 4 . In Eq. JT3JI the relative strength 
of the singularity also depends on the distance between 
the QCD critical end point, where S c — 0, and the 
s q = sh line. If the CEP is located near the sq = sh 
line, less singularity is realized. However, in the follow- 
ing calculations, we assume that the CEP is on the phase 
transition line in the bag plus excluded volume model, 
where the phase transition is always of strong first order 
as shown later. Thus, sq is substantially larger than sh 
at the CEP and the singularity is not suppressed. 

Once we construct the entropy density, we can calcu- 
late the other thermodynamical quantities such as the 
baryon number density, pressure, and energy density. 
The baryon number density ns is given by 



n B (T, fi B ) 



dP 
dj-iB 
T ds(T',fi B ) 
dfi B 



dT' + n B (0, f i B ).(16) 



P(T,fi B )= [ s(T\fi B )dT' +P(0,mb). 
Jo 



(18) 



Finally, from the thermodynamical relation the energy 
density E is obtained as 



Ts - P - hbtib- 



(19) 



Here we make a comment on the choice of the param- 
eters AT cr it, A/iscrit, and D. The linear mapping from 
(r, h) to (T, /xb), strictly speaking, holds only in the prox- 
imity of the CEP, where the thermodynamical quantities 
comply with the universality. We use the linear map- 
ping throughout the critical region. In order to avoid 
artifacts due to this simplification, AT cl .j t , A/i Bcr i t , and 
D need to be appropriately chosen so that thermody- 
namical constraints are not violated. In particular, the 
following thermodynamical inequalities [23, 124I |25| need 
to be observed, 



U1 / V,N B 



> 0, 



T,N n 



dP 
dV 

^) > 



(20) 



where S is the entropy, Nb is the baryon number, and V 
is the volume of the system. 

Figure |3 indicates the entropy density and the baryon 
number density as a function of T and fiB- We can see 
that in the high chemical potential region the order of 
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the phase transition is first and that in the low chemical 
potential region a smooth phase transition occurs, which 
reflects the existence of the CEP. In this calculation the 
CEP is located at (T,n B ) = (154.7 MeV, 367.8 MeV). 
Since the universality does not provide the information 
on the location of the CEP, it is treated as a parame- 
ter throughout this paper as well as the values of AT cr ;t, 
A/xscrit, and D. Note, however, that the local singular 
behavior of thermodynamical quantities around the CEP 
is fixed by the universality hypothesis and that accord- 
ingly the local features such as the focusing, which we 
discuss below, are not affected by this ambiguity. 




He [MeV] 



FIG. 3: Entropy density (left) and the baryon number 
density (right) as a function of T and \ib- (Te,/ibe) = 
(154.7 MeV, 367.8 MeV). The values of (AT crit , A^ Bc rit, D) 
are the same as in Fig. 0J (left). 



Two comments are in order here. First, as we ex- 
plained, the universality does not tell us about the sizes 
of AT cr j t , A/iB cr it, D, and so forth. If D is larger, the size 
of the critical region becomes smaller and the behavior 
of the isentropic trajectories approaches that in Fig. |SJ 
which can be seen from Eqs. (|12|l and 1)13(1. Namely in 
the case of larger D, the contour lines of S c (for exam- 
ple, S c = —1, —0.5,0.5 and 1 in Fig. are closer to the 
S c = line, i.e., the tangential line to the QCD phase 
boundary at the CEP. Then, in the limit of D — » oo, 
all contour lines shrink to one line, which is identical to 
the QCD phase boundary as long as the curvature of 
the phase boundary can be ignored. This can be seen 
from Eq. ifT^I . Therefore, the bag plus excluded volume 
model can be considered as an extreme case of the equa- 
tion of state with the CEP. Second, the focusing on the 
right hand side of the CEP was first discussed in Ref. 
|16|. Our numerical result indicates that the CEP at- 
tracts isentropic trajectories not only on the right hand 
side (first order side) of the CEP but also on the left 
hand side (crossover side) of the CEP. Thus, if the criti- 
cal region is large enough, the fine-tuning of the collision 
energy is not necessary to hit the vicinity of the CEP 
not only on the low energy side [l6l | but also on the high 
energy side. 



First we investigate the behavior of the isentropic tra- 
jectories, i.e., contour lines of ub/s. When entropy pro- 
duction can be ignored, the entropy and baryon number 
are conserved in each volume element and, therefore, the 
temperature and chemical potential in a given volume 
element change along the contour lines specified by the 
initial condition. 

Figure 0| shows the isentropic trajectories in the T-\xb 
plane. The values of (AT crit , A/i Bcri t, D) for the Fig.0] 
(left) and Fig. H (right) are (100 MeV, 200 MeV, 0.15) 
and (100 MeV, 200 MeV, 0.4), respectively. In both 
cases, the trajectories are focused to the CEP. This 'fo- 
cusing' effect is more clearly observed in the case of 
(T e ,mbe) = (143.7 MeV, 652 MeV) (Fig. H (right)). 
Thus, the CEP acts as an attractor of isentropic trajec- 
tories. Figure [U shows isentropic trajectories in the bag 
plus excluded volume model, which is currently employed 
in most of hydrodynamical calculations |26(. The order 
of the phase transition is always first in this case. There 
is no focusing effect on the isentropic trajectories in this 
case. Instead, the trajectories are just shifted to the left 
on the phase transition line. This implies that the hy- 
drodynamical evolution in the case with the CEP is very 
different from the one in the case with the equation of 
state given by the bag plus excluded volume model. Be- 
cause of this attractor character of the CEP, it is not 
needed to fine-tune the collision energy to make the sys- 
tem pass near the CEP, which is pointed out in Ref. [R| 
for the first time. In other words, it is expected that 
the effect of the CEP, if any, changes only slowly as the 
collision energy is changed. 
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FIG. 4: Isentropic trajectories in the cases with the CEP. The 
CEP is located at (T e ,mbe) = (154.7 MeV, 367.8 MeV) 
(left) and (T e ,/xb E ) = (143.7 MeV, 652.0 MeV) (right). 
The values of ub/s on the trajectories are 0.01, 0.02, 0.03, 
and 0.04 (left) and 0.01, 0.02, 0.03, 0.04, and 0.05 (right) 
from left to right. L/L to tai(T, fis) is shown for some points 
on two trajectories in the left figure. 

Next we argue the change of the square of the sound 
velocity c^, 



dP 
dE 



(21) 



n B /s 



along the isentropic trajectories, which gives us the in- 
formation on how the equation of state changes in the 
course of the hydrodynamical evolution. Figure |H1 indi- 
cates the sound velocity as a function of L/L to tai, where 
L is the path length to a point along the isentropic tra- 
jectory with a given ub/s from a reference point on the 
same isentropic trajectory on the T-fis plane, and Ltotai 
is the one to another reference point along the trajectory. 
For example, on the insentropic trajectory ng/s = 0.01 
in Fig. |U (left), L/L tota \ — 0, which is the starting point 
of ub/s trajectory, corresponds to the point (T,/is) = 
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FIG. 5: Isentropic trajectories (solid lines) in the bag plus 
excluded volume model. The values of ns/s on the trajecto- 
ries are 0.01, 0.02, 0.03, 0.04, and 0.05 from left to right. The 
dashed line stands for the phase boundary. The order of the 
phase transition is always first in this case. 



(200 MeV, 145 MeV) and L/L tota x = 0.29 corresponds 
to the point (T,fj, B ) = (155 MeV, 361 MeV). At about 
L/L tota i = 0.25 ~ 0.35 the hadronization process occurs 
and at about L/Xtotai = 0.35 ~ 0.5 kinetic freezeout 
takes place. We can see clear difference between the be- 
havior of the sound velocities on ns/s = 0.01 and 0.03 
trajectories. In the case of ub/s = 0.01, the sound ve- 
locity changes smoothly along the trajectory, which re- 
flects the occurrence of a smooth phase transition at this 
ns/s value, and it takes its minimum at L/L to tai = 0.29. 
On the other hand, the sound velocity at ns/s = 0.03 
changes suddenly at L/L to ta\ = 0.34, due to the first 
order phase transition. The hydrodynamical expansion 
along various tib/s paths differs in how the effect of the 
phase transition appears because of the existence of the 
CEP. In a real collision system, ns/s changes from posi- 
tion to position. As a result, the time evolution of such 
a system is described as a superposition of trajectories 
with different ng/s's. Thus, the system may not expand 
uniformly and this effect may appear in physical observ- 
ables related to expansion such as the collective flow and 
Hanbury Brown- Twiss (HBT) radii if the collision pa- 
rameters are appropriate, i.e., isentropic trajectories in 
the system pass through and/or near the CEP. 
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III. SLOWING OUT OF EQUILIBRIUM 

The equilibrium correlation length and fluctuation be- 
come large near the CEP and they diverge at the CEP. At 
the same time, the typical time scale becomes long near 
the CEP, which makes the system take long time to reach 
the thermal equilibrium near the CEP. Therefore, in the 
time evolution in ultrarelativistic heavy ion collisions, 
near the CEP slowing-out-of-equilibrium occurs and non- 
equilibrium dynamics has to be taken into account. We 
assume that the thermal equilibrium is established soon 
after collisions and that the quark-gluon plasma follows 
Bjorken's scaling solution |27| for hydrodynamical evolu- 
tion. The initial temperature and proper time are set to 
200 MeV and 1 fm/c, respectively 2 . Before discussing 
the correlation length in out-of-equilibrium time evolu- 
tion, we calculate the equilibrium correlation length near 
the CEP in the r-h plane. Using Widom's scaling law 
|28j . the equilibrium correlation length £ cq is given by 



a 2 M -2u/0 g 



IMIV/3 



(22) 



where v — 0.63 is a critical exponent of the three dimen- 
sional Ising model |2l|. a is a constant with the dimen- 
sion of length and fixed to 1 fm. There is not so large 
ambiguity in the determination of a, since £ e q is of the 
order of 1 fm at T ~ 200 MeV The function g(x) is given 
by the e expansion to order e 2 psj . 



g(x) = g e (x) 

= &- 2y z{\-— [(5 + 61n3)z-6(l + z)lnzl 
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In 3 



_ fl ln[a(l-z)] 



137 

dx ~ 



3 J 

(23) 

-2.344, and 



where z = . , 2 jQ , I = L — r . 1 

l+K/S' JO 1— x{l— x) 

e = 4 — d with d being the dimension of the space, 3. At 
large x, i.e., around the crossover and CEP, g e (x) cannot 
be used and the asymptotic form given in Ref. [29j, 



g(x) = gi avse (x) = hrr- - ( 24 ) 



FIG. 6: Square of the sound velocity c s as a function of 
L/ Ltotal along ns/s = 0.01 (left) and ns/s = 0.03 (right) 
isentropic trajectories. The CEP is located at (Te,^be) = 
(154.7 MeV, 367.8 MeV) in both cases. The parameters are 
the same as in Figs. 2, 3, and 4 (left). 



2 In Bjorken's one-dimensional hydrodynamical evolution, the 
non-equilibrium correlation length is not sensitive to the choice 
of the initial conditions. Instead, as we discuss later, the reduc- 
tion of the time scale due to the transverse expansion is expected 
to have more effect on it. 
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should be used. We smoothly connect g e (x) and siarge(^) 
around x = 5, i.e., in most part of the critical region g{x) 
is given by Eq. jSBjl - 

Figure shows the correlation length in equilibrium 
and the isentropic trajectories at some ng/s's in the r-h 
plane. In Fig. the CEP is located at the origin, and 
the first order phase transition occurs at r < and the 
crossover phase transition takes place at r > 0. The 
correlation length in equilibrium £ eq is divergent at the 
origin and the region where £ cq is large spreads out in 
the smooth phase transition region, r > 0. In Fig. [7| 
solid lines stand for the isentropic trajectories in the 
case of (T E ,fi BE ) = (154.7 MeV, 367.8 MeV) (left) and 
(T e ,Hbe) = (143.7 MeV, 652.0 MeV) (right). Note 
that this r-h plane is actually placed in the T-fi B plane 
with the tilt given by the slope of the tangential line of 
the first order phase transition line at the CEP as shown 
in Fig. The value of n B / s is determined not only by 
r and h but also by the non-critical component of the 
equation of state. This is the reason why the behavior 
of the trajectories relative to the CEP changes according 
to the location of the CEP in the T-/i B plane. Thus, 
the strength of focusing effect of the CEP is determined 
also by the non-singular part. On the other hand, we 
use the universality hypothesis to calculate £ e q and, as a 
result, it is a function of only r and h. The trajectories 
are attracted to the CEP, which makes the length of each 
trajectory in the critical region, where £ oq is large, long. 
This behavior is very different from the assumption in 
usual schematic analyses, for example, Fig. 1 in Ref. [2£j 
by Berdnikov and Rajagopal. The maximum value of £ e q 
depends on the value of n E /s. 



\aoi 
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1st order mU.o- 




m >2.s 







FIG. 7: Correlation length in equilibrium and the tib/s 
contour lines in the r-h plane. The solid lines are trajec- 
tories with fixed ns/s values 0.008 and 0.01 in the case of 
{Te,hbe) = (154.7 MeV, 367.8 MeV) (left) and 0.02 and 
0.035 in the case of {T e ,^be) = (143.7 MeV, 652.0 MeV) 
(right). 

Next we calculate the non-equilibrium correlation 
length £(t) as a function of the proper time r by using 
the rate equation given in Ref. |29j , 



— m ff (r) = -r[nv(T)] ( m CT (r) - 



(25) 



r[m cr (r)] is given as 

r Kr(r)] = ^(m CT (T)£ ) 



(26) 



where 1/to (T (t) = £(r) and r[rn cr (T)] is the parameter 
which represents the rate of slowing-out-of-equilibrium. 
From the theory of dynamical critical phenomena, 



where we take £o = £ e q(T = 200 MeV), and z is a univer- 
sal exponent of Model H in |3(j, i.e., z ~ 3 p9ll3l| . A is a 
dimensionless non-universal parameter and we use A = 1 
as in Ref. |2^, though there is considerable uncertainty 
in the parameter A. However, fortunately, the relation 
between £ e q and £ does not depend on the parameter 
choice of A so much at the kinetic freezeout point. The 
non-equilibrium correlation length approaches the equi- 
librium correlation length as A increases. At A — 100, £ 
is almost equal to £ cq . On the other hand, at small A the 
effect of slowing-out-of-equilibrium appears strongly. At 
the same time, however, the maximum value of £ becomes 
small as A decreases. As a result, the absolute value of 
the difference between £ cq and £ at the kinetic freezeout 
temperature remains almost unchanged for a wide range 
of A. Indeed, we have checked that the following discus- 
sion about non-equilibrium correlation length holds for 
0-1 < A < 10. 

Figure [S] shows the correlation length as a function 
of L/Ltotai- The dashed and solid lines stand for the 
correlation lengths in equilibrium at n E /s = 0.008 and 
n B /s — 0.01, respectively. The parameters are the same 
as for Fig. 01 (left) . The thin and thick lines are £ e q and £, 
respectively. The maximum value of £ oq along the former 
trajectory is larger than that along the latter, because the 
former approaches the CEP more closely than the latter. 
The non-equilibrium correlation length £ is smaller than 
£ oq at the beginning. Then, £ becomes larger than £ e q 
later. These are both due to the critical slowing down 
around the CEP, as pointed out in Ref. 29] . However, 
the difference becomes small by the time the system gets 
to the kinetic freezeout point. If the transverse expansion 
is taken into account, the time scale in the hadron phase 
becomes much shorter, but £ e q is already small in the 
hadron phase and the difference is expected to remain 
small. Figure |21 shows the cube of the correlation length 
as a function of the temperature for n E /s = 0.008, 0.01, 
and 0.015. The parameters are the same as for Fig. [S] 
The fluctuation of the sigma field scales approximately 
as oc £ 3 . The singular part of the fluctuation of intensive 
thermal quantities of pions near the CEP is proportional 
to £ 2 ^(|. The general relation between £ and observed 
fluctuations still needs to be clarified. These relations, 
however, suggest that there is room for the observation 
of some enhancement in fluctuations when the system 
passes the vicinity of the CEP even if the enhancement 
of £ is small. 



IV. CONSEQUENCES OF THE CRITICAL END 
POINT IN EXPERIMENTS 

Many experimental analyses have been carried out in 
search of the evidence of the CEP in the QCD phase 
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FIG. 8: Equilibrium correlation length £ cq (thin lines) and 
non-equilibrium correlation length f; (thick lines) on the isen- 
tropic trajectories with ub/s = 0.008 and 0.01, together with 
r as functions of L/L tota ,i (inlet), z = 3 was used in the 
calculation. 
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FIG. 9: Cube of the equilibrium correlation length (thin 
lines) as a function of the temperature and the cube of the 
non-equilibrium correlation length (thick lines) on the isen- 
tropic trajectories with ub/s = 0.008, 0.01, and 0.015. The 
parameters are the same as for Fig. |S] The isentropic lines 
with ub/s = 0.008, 0.01, and 0.015 pass left of, almost 
through, and right of the CEP, respectively. 

diagram. Since the correlation length and fluctuation di- 
verge at the CEP in thermal equilibrium, it has been 
expected that some enhancement of the fluctuation, for 
instance, is observed if the collision energy of nuclei is 
properly adjusted so that the system goes right through 
the CEP. Generally, the higher the collision energy is, 
the smaller chemical potential region is explored in the 
T-/ib plane. Thus, it has also been expected that the 
observables related to the critical behavior around the 
CEP such as fluctuations show non-monotonic behavior 
as a function of the collision energy. We discuss the ob- 
servability of such behavior in ultrarelativistic heavy ion 
collisions on the basis of our findings in the previous sec- 
tions. 

It has been naively expected that the collision energy 
needs to be carefully adjusted so that the system goes 
right through the CEP and its existence can be con- 
firmed. In section II, we have shown that the CEP acts 




100 200 300 400 500 600 
Hb [MeV] 



FIG. 10: Isentropic trajectories and entropy density contour 
lines in the case with the CEP. The parameters are the same 
as for Fig. [3] and the left figure of Fig. 2] The solid lines are 
isentropic trajectories (n B /s = 0.01, 0.015, 0.02, 0.025, 0.03, 
0.035, and 0.04 from left to right) as in Fig. The dashed 
lines are entropy density contour lines, which are shown at 
every 1 fm -3 . 



as an attractor of the isentropic trajectories. Thus, if the 
size of the critical domain is large enough, as stressed in 
Ref. 0| for the right hand side of the CEP, it is not neces- 
sary to fine-tune the collision energy to make the system 
approach the CEP closely enough. However, such phys- 
ical quantities that show the critical behavior near the 
CEP are hadronic observables, and are subject to final 
state interactions in the hadron phase. It has been often 
argued that if the system passes near the CEP, kinetic 
freezeout takes place near the CEP [l^. It is based on 
the expectation that, in such a case, the system stays long 
near the CEP, where the phase transition is second order, 
and that when the system starts to leave away from the 
CEP, the entropy density or particle density is already 
small enough so that kinetic freezeout takes place. To 
the contrary, Fig. [3] does not show a sharp drop of the 
entropy density near the CEP. To see this more clearly, 
we show the contour plot of the entropy density in Fig. 1101 
for the same parameters as in Fig. 01 The entropy density 
gives a semi-quantitative measure of the whereabouts of 
kinetic freezeout. Alternatively, we could use the energy 
density as a measure as well. Figure ITU1 shows that the 
contour lines of the entropy density are not focused near 
the CEP in contrast to the isentropic trajectories. Thus, 
it is not likely that kinetic freezeout takes place near the 
CEP. The reason can be traced back to the fact that at 
the CEP the entropy density is continuous and bound, 
although its derivative along the h axis diverges. If the 
system does not freeze out near the CEP, the dilution 
of the critical behaviors in the hadron phase needs to be 
carefully considered. 

Empirically it has been known that chemical freezeout 
takes place just below the theoretically expected phase 
transition line. Since the isentropic trajectories are fo- 
cused around the CEP, it would be tempting to expect 
that the chemical freezeout points are focused as the col- 
lision energy is varied. However, so far, only the free 
resonance gas model has been used in the analysis of 
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the chemical freezeout point. If the critical region is so 
large that the focusing is realized in a large region, the 
effect of the interactions on hadrons can by no means 
be ignored. Then, the hadrons at chemical freezeout are 
at most quasi-particles which are quantum-mechanically 
different from the hadrons in the vacuum |32j | and it is 
necessary to consider the time evolution of the quasi- 
particle states to relate the chemical freezeout points and 
the observed particle number ratios. This is beyond the 
scope of this paper and we would like to leave this as a 
future problem. 

The entropy density in Fig. [21 shows the following fea- 
ture. On the left hand side of the CEP, the entropy den- 
sity changes smoothly, reflecting the fact that the phase 
transition here is crossover. The closer to the T axis, 
the less rapid the change is. On the other hand, on the 
right hand side of the CEP, the phase transition is of first 
order, and the jump in the entropy density increases as 
the chemical potential increases. As a result, i) the sep- 
aration of the contour lines of the entropy density and 
energy density gets larger as the distance from the CEP 
becomes larger on the left hand side of the CEP, ii) the 
contour lines below the crossover are flatter than those 
in the hadron phase in the case with a first order phase 
transition (See Figs. l3l and 110(1 . This feature is not due 
to the universality. Instead, it is a general consequence 
of the existence of the CEP. In fact, as shown in Fig. 1111 
this feature is not observed in the case of the bag plus ex- 
cluded volume model. Using this feature of the entropy 
density, we discuss the behavior of kinetic freezeout as a 
function of the collision energy. Kinetic freezeout takes 
place when the mean (elastic) collision rate and the ex- 
pansion rate of the system become comparable. The ex- 
pansion rate in central collisions at RHIC is known to 
be considerably larger than that at SPS [33j. The elastic 
collision rate can be estimated from the elastic cross sec- 
tions for pions and baryons and their densities. Due to 
the chiral symmetry and the effect of the A, the elastic 
ir-N cross section is considerably larger than that of tt-it. 
On the other hand, in Pb+Pb collisions at 158 AGeV at 
SPS, (p + v)/(n+ +7T-) - 0.09 and it is - 0.09 at 

RHIC 36] . Thus, at the same particle number density, 
the collision rate at SPS is almost the same as that at 
RHIC. From these considerations, it would be natural if 
the kinetic freezeout temperature in central collisions at 
RHIC is noticeably higher than that at SPS, while the ki- 
netic freezeout temperature at RHIC is actually slightly 
lower than that at SPS |33|. Two comments are in or- 
der here. First, the kinetic freezeout we are conside ring 
here is that of pions and kaons. The analysis in Ref. |33| 
does not take into account the effect of resonance de- 
cays. When resonance decays are taken into account, the 
kinetic freezeout temperature becomes even lower |37l . 
Second, more precisely speaking, kinetic freezeout takes 
place when the mean collision rate and the product of 
the expansion rate and the typical scale of the system 
become comparable. However, according to the results 
of the HBT interferometry, the typical size at RHIC is 




MMeV] 

FIG. 11: Isentropic trajectories and entropy density contour 
lines in the case without the CEP. The parameters are the 
same as in Fig. The solid lines are isentropic trajectories 
(n B /s = 0.01, 0.02, 0.03, 0.04, and 0.05 from left to right) as 
in Fig. |K| The dashed lines are the entropy density contour 
lines. The contour lines are shown at every 0.1 (1) fm -3 in the 
hadron (QGP) phases because the increase rate of the entropy 
density in the QGP phase is much larger than that in the 
hadron phase due to the strong first order phase transition. 



not increased as the expansion rate compared to SPS. 

This can be naturally understood if the separation of 
the entropy density (<~ particle density) or energy den- 
sity contour lines is larger at kinetic freezeout at RHIC 
than at SPS. This implies that the CEP exists in the 
QCD phase diagram and, furthermore, that the RHIC 
isentropic trajectory passes on the left hand side of the 
CEP. To confirm this tendency, it is necessary to obtain 
the kinetic freezeout temperatures at several collision en- 
ergies between the SPS and current RHIC energies. 

In the previous section, we have shown that the differ- 
ence between the thermal correlation length £ cq and the 
actual correlation length £ is small at kinetic freezeout 
as well as £ itself. As explained in Ref. ;29j, precisely 
speaking, Eq. I|26|) is valid only near the CEP. Neverthe- 
less, the physical reason why £ eq and £ are so close to 
each other at kinetic freezeout, is general and holds also 
in this case; final state interactions tend to wash out non- 
equilibrium effects. Since the kinetic freezeout points at 
the Super Proton Synchrotron (SPS) at CERN and the 
Relativistic Heavy Ion Collider (RHIC) at Brookhaven 
National Laboratory are far below the phase transition 
line, the effect of final state interactions is expected to 
dominate. 

Similar arguments hold also for fluctuations except for 
the cube effect we argued in the previous section. The 
event-by-event fluctuations of the mean transverse mo- 
mentum in Pb + Au collisions at 40, 80, and 158 AGeV/c 
were measured by the CERES Collaboration |38|. They 
are slightly smaller than those at RHIC, and unusually 
large fluctuation or non-monotonic behavior which may 
suggest the existence of the CEP has not been observed. 
However, in order to understand experimental observa- 
tions, we have to take account of the CEP character as 
an attractor of isentropic trajectories in the T-/ib plane. 
Due to the focusing effect of the CEP, the isentropic tra- 
jectories of various initial collision energies are gathered 
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to the CEP, if the critical region around the CEP is large 
enough. As a result, when each isentropic trajectory 
passes near the CEP, similar correlation lengths, fluctua- 
tions, and so on are induced, and they do not show strong 
non-monotonic behavior as a function of the collision en- 
ergy. In other words, the absence of the non-monotonic 
behavior, which is shown by the CERES Collaboration, 
does not necessarily imply the non-existence of the CEP 
in the region probed by SPS and RHIC. More precise 
measurement and analyses which take account of the fo- 
cusing and slowing-out-of-equilibrium may reveal the ex- 
istence of the CEP. 

Recently, the fluctuation of (K + + K~)/(tt + + n~) 
and (p + p)/(tt + + n~) was measured and it was found 
that the collision energy dependence of the fluctuation 
of (K + + K~)/(ir + + 7r~) is not described by a cascade 
model 34]. While this result is interesting, we have to 
be careful in associating this result with the critical phe- 
nomena around the CEP. Since neither the numerator 
nor denominator of (K + + K~)/ (it + +7r _ ) is a conserved 
quantity in the strong interaction, substantial modifica- 
tion is expected for both in the hadron phase as we have 
demonstrated for the correlation length. Furthermore, as 
we discussed for the particle ratios, quasi-particle states 
near the CEP are quantum-mechanically different from 
the ones in the vacuum, i.e., observed pions, kaons, pro- 
tons, and so on. To relate the fluctuation of the particle 
ratios to the critical phenomena, it is necessary to know 
not only the time evolution of the quasi-particle states 
but also the relation between the quasi-particle quantum 
states in medium and the particle states in the vacuum 

m 

Up to now, we have been considering the isentropic 
trajectories in the T-/ib plane. Actually, the viscosity 
is known to diverge at the CEP and the entropy is gen- 
erated if the system passes near the CEP. However, the 
singularity in the viscosity is known to be very weak [3l| , 



Scq 
£cq,0 



(27) 



where 77 is the shear viscosity, r/o is the shear viscosity 
at the point where the equilibrium correlation length is 
£ cq ,o, and e is 4 — d as before. In Fig. El v/vo in the 
r-h plane is shown. The small correction to the criti- 
cal exponent 1/19 is set to for simplicity. Because of 
the small critical exponent, the diverging feature of the 
viscosity is hardly seen. The effect of the entropy pro- 
duction around the CEP is thus expected to be small, 
and it will not show easily-recognizable non-monotonic 
changes in observables as the collision energy is changed. 



V. SUMMARY 

In this paper, we discussed the hydrodynamical expan- 
sion near the QCD critical end point and its consequences 




FIG. 12: Shear viscosity rj/rjo near the CEP. For the param- 
eters used in the calculation, see the text. 



in experimental observables. First we constructed realis- 
tic equations of state with the CEP on the basis of the 
universality hypothesis and discussed the hydrodynami- 
cal expansion near the CEP. The behavior of isentropic 
trajectories near the QCD critical end point is clearly dif- 
ferent from that in the bag plus excluded volume model, 
which is usually used in hydrodynamical models. We 
found that the CEP acts as an attractor of isentropic 
trajectories, Ub/s = const., in the T-fiB plane, not only 
on the right hand side of the CEP but also on the left 
hand side of the CEP. Because of the focusing of the isen- 
tropic trajectories, the path of the system in the course 
of the time evolution in the T-^b plane is remarkably dif- 
ferent from that in the usual bag plus excluded volume 
model. This will be reflected in the final state observ- 
ables that are sensitive to the equation of state, such as 
the anisotropic flows (directed, elliptic, and higher order 
Fourier coefficients) [23, and HBT radii. For this pur- 
pose, full 3-D hydrodynamical calculations I4ll[42l with 
realistic equations of state are indispensable |43j |. 

Next we argued the critical slowing down near the 
QCD critical end point. We compared the equilibrium 
and non-equilibrium correlation lengths along two isen- 
tropic trajectories which pass near the CEP. We found 
that the difference between them at kinetic freezeout is 
very small. 

Furthermore, we considered the experimental observ- 
ability of the consequences of the CEP. Fluctuations and 
viscosities diverge at the CEP, but fluctuation which is 
induced by the CEP will fade by kinetic freezeout for the 
same reason as for the correlation length. This explains 
the non-observation of sharp non-monotonic behaviors in 
fluctuations. However, the singular part of the fluctua- 
tions is approximately proportional to the cube of the 
correlation length and there is a possibility that more 
precise measurements in the future may reveal the ex- 
istence of the CEP. The singularity in the viscosity is 
too weak to affect observables. Concerning the chemical 
freezeout process, we need to consider the effect of the 
interactions on hadrons for detailed discussion on parti- 
cle number ratios, when it takes place near the CEP. We 
found that, at least now, the most promising way is to 
trace the behavior of the collision energy dependence of 
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the kinetic freezeout temperature. The fact that the ki- 
netic freezeout temperature at RHIC is lower than that 
at SPS suggests the existence of the CEP in the T-fis 
plane and that the RHIC path goes left of the CEP. 

Although our discussion in this paper does not depend 
on the parameter choice in Eqs. (|12(l and (|13(l qualita- 
tively, our assumption that the critical region near the 
CEP spreads out a large area may look too optimistic. 
However, according to recent lattice calculations at zero 
baryon chemical potential, the ipip susceptibility shows a 
very broad peak as a function of the temperature. Fur- 
thermore, the width tends to even increase as the lattice 
spacing is decreased |44| . This clearly indicates that the 
free-gas description cannot be used at 100 < T < 250 
MeV at zero baryon chemical potential. 



lowing relations [24l l25| , 



SS = 0, (Al) 
S 2 S < 0, (A2) 



are satisfied. Here the independent variables are the en- 
ergy E, volume V, and particle number N. Suppose 
that the total system whose volume is V is composed of 
two parts with the volumes V\ and Vi. Suppose that in 
V% {V 2 ) the variation of the energy, volume, and parti- 
cle number are SEi (SE 2 ), 5V\ (SV 2 ), and SNi (5N 2 ), 
respectively. These variations satisfy 
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APPENDIX A: STABILITY CONDITION 

In an isolated system, the entropy S takes the max- 
imum value in the stable equilibrium state and the fol- 



SE 1 + SE 2 = 0, 
5V1 + 8V2 = 0, 

5Ni + 5N 2 = 0. (A3) 



By expanding SS with regard to E, V , and N in Eq. 
dSU and substituting (ff) yjv = (w)e,n = £> 

(§w)ev = ~T' and E( i- ( E3 int0 E i- E2J> we ob_ 

tain T\ = T 2l P\ — P 2l and [i\ = fi 2 , i.e., both parts 
of the system have the same temperature, pressure, and 
chemical potential in equilibrium. 

Next we evaluate the second variation 5 2 S in Eq. I|A2(I . 



S 2 S 



2 S 1 



S 2 E 1 



d 2 s 1 



+2 



5 2 S 1 



SExSVx 
+ (1-2). 



8 2 V X 
5 2 S, 



d 2 s 1 



5^ ■ SVi + 2 



S 2 ^ 



S 2 S 1 
SEiSNi 



SEi ■ 5Ni 



Vi 



(A4) 



By calculating S 2 S using T5S = 5E + PSV — fiSN and By definition, 
comparing it with (|A4p . we obtain 

T5 2 S = -STrfSi + SPiSVi-SmSNt + (1 -> 2). (A5) 
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SSi = 
SPi = 
Sjj,i = 



asi 

dP x 
dTi 



^1 <m 

^ ) ^ 



5Vi 



dSi 

w) sv ' 

dV lJT uNl 



V 9iVi 



0Pi 
dm 



8N lt 



T u Vx 



TuVt 



SNi. 



(A6) 



T u Vi 



Substituting Eq. (|A6|) and the following relations, and using Eq. (| A3|> . we obtain the quadratic form 

\dT)V,N T v ' \dV)T,N \dT)V,N> \ dV J T N 



5V 



T.N 



(M.) 

\dN )T,V 

\dV JT,N 



SN 



dN . 



( dP\ 2 
KdNlT.y 

TIE] 

T,V \dV>T,N 



(A7) 



where we have suppressed the suffixes on the right hand We note that (|A10|I is equivalent to the following inequal- 
side. Finally from Eqs. (|A2|) and l|A7[l . we obtain the ity |24| . 
stability conditions, 

C v > 0, (A8) 



T,N 

2 



^ >-jf^>0. (A10) 



dN 
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